{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this notebook we will compare bnlearn's hillclimbing algorithm and tabu algorithm against Gobnilp. We'll do this with default settings and BIC scoring. First let's use bnlearn's hill-climbing algorithm to learn from its built-in alarm dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n",
      "Attaching package: ‘bnlearn’\n",
      "\n",
      "The following object is masked from ‘package:stats’:\n",
      "\n",
      "    sigma\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "  Bayesian network learned via Score-based methods\n",
       "\n",
       "  model:\n",
       "   [HIST][HRBP][PAP][FIO2][ANES][ERCA][LVF|HIST][PMB|PAP][ERLO|HRBP][PCWP|LVF]\n",
       "   [HR|HRBP:ERLO][HREK|HR:ERCA][HRSA|HR:ERCA][LVV|PCWP:LVF][CCHL|HR][CVP|LVV]\n",
       "   [MINV|CCHL][STKV|LVF:LVV][CO|STKV:HR][HYP|LVV:STKV][VALV|MINV][INT|MINV:VALV]\n",
       "   [PVS|FIO2:VALV][ACO2|CCHL:VALV][PRSS|INT:VALV][SHNT|PMB:INT]\n",
       "   [VLNG|MINV:INT:VALV][SAO2|SHNT:PVS][ECO2|ACO2:VLNG][KINK|PRSS:VLNG]\n",
       "   [VTUB|PRSS:MINV:INT][TPR|SAO2:CCHL][DISC|VTUB][BP|TPR:CO][APL|TPR]\n",
       "   [VMCH|DISC:VTUB][MVS|VMCH]\n",
       "  nodes:                                 37 \n",
       "  arcs:                                  53 \n",
       "    undirected arcs:                     0 \n",
       "    directed arcs:                       53 \n",
       "  average markov blanket size:           3.46 \n",
       "  average neighbourhood size:            2.86 \n",
       "  average branching factor:              1.43 \n",
       "\n",
       "  learning algorithm:                    Hill-Climbing \n",
       "  score:                                 BIC (disc.) \n",
       "  penalization coefficient:              4.951744 \n",
       "  tests used in the learning procedure:  2718 \n",
       "  optimized:                             TRUE \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "-220761.687712649"
      ],
      "text/latex": [
       "-220761.687712649"
      ],
      "text/markdown": [
       "-220761.687712649"
      ],
      "text/plain": [
       "[1] -220761.7"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "library(bnlearn)\n",
    "data(alarm)\n",
    "alarmhc.bn <- hc(alarm)\n",
    "alarmhc.bn\n",
    "score(alarmhc.bn,alarm,type=\"bic\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "OK, so hill-climbing has quickly found a BN with BIC score of -220761.687712649. Let's see how tabu search does."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "  Bayesian network learned via Score-based methods\n",
       "\n",
       "  model:\n",
       "   [PAP][FIO2][LVF][ANES][ERLO][HR][ERCA][HIST|LVF][HRBP|ERLO:HR][HREK|HR:ERCA]\n",
       "   [HRSA|HR:ERCA][PMB|PAP][LVV|LVF][CCHL|HR][CVP|LVV][PCWP|LVV][MINV|CCHL]\n",
       "   [STKV|LVF:LVV][CO|STKV:HR][HYP|LVV:STKV][VALV|MINV][INT|MINV:VALV]\n",
       "   [PVS|FIO2:VALV][ACO2|CCHL:VALV][PRSS|INT:VALV][SHNT|PMB:INT]\n",
       "   [VLNG|MINV:INT:VALV][SAO2|SHNT:PVS][ECO2|ACO2:VLNG][KINK|PRSS:VLNG]\n",
       "   [VTUB|PRSS:MINV:INT][TPR|SAO2:CCHL][DISC|VTUB][BP|TPR:CO][APL|TPR]\n",
       "   [VMCH|DISC:VTUB][MVS|VMCH]\n",
       "  nodes:                                 37 \n",
       "  arcs:                                  51 \n",
       "    undirected arcs:                     0 \n",
       "    directed arcs:                       51 \n",
       "  average markov blanket size:           3.41 \n",
       "  average neighbourhood size:            2.76 \n",
       "  average branching factor:              1.38 \n",
       "\n",
       "  learning algorithm:                    Tabu Search \n",
       "  score:                                 BIC (disc.) \n",
       "  penalization coefficient:              4.951744 \n",
       "  tests used in the learning procedure:  3759 \n",
       "  optimized:                             TRUE \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "-220727.331509064"
      ],
      "text/latex": [
       "-220727.331509064"
      ],
      "text/markdown": [
       "-220727.331509064"
      ],
      "text/plain": [
       "[1] -220727.3"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "alarmtabu.bn <- tabu(alarm)\n",
    "alarmtabu.bn\n",
    "score(alarmtabu.bn,alarm,type=\"bic\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tabu search is also very quick and finds a slightly better network with a score of -220727.331509064.\n",
    "\n",
    "Now let's see how Gobnilp does on the same data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "**********\n",
       "BN has score -218632.31567964845\n",
       "**********\n",
       "ACO2<-CCHL,VALV -7365.911141511497\n",
       "VALV<-INT,VLNG -3656.8027442121847\n",
       "CCHL<-SAO2,TPR -5079.249250888837\n",
       "ANES<- -6708.582939203498\n",
       "APL<- -1134.1626027276418\n",
       "BP<-CO,TPR -10757.751590048356\n",
       "TPR<-APL -21601.664138686014\n",
       "CO<-HR,STKV -5340.866593873159\n",
       "SAO2<-PVS,SHNT -2469.044793071091\n",
       "HR<-CCHL -7361.101683937812\n",
       "STKV<-HYP,LVF -9121.2931288759\n",
       "CVP<-LVV -6151.495470079891\n",
       "LVV<-HYP,LVF -7523.033824999751\n",
       "DISC<- -6447.082832428477\n",
       "ECO2<-ACO2,VLNG -3540.04821713251\n",
       "VLNG<-INT,VTUB -7559.536030079053\n",
       "ERCA<- -6480.204445349012\n",
       "ERLO<- -3981.1433847172025\n",
       "FIO2<- -4007.581975260592\n",
       "HIST<-LVF -1488.3560086436557\n",
       "LVF<- -3945.7594222796324\n",
       "HRBP<-ERLO,HR -2822.433256557475\n",
       "HREK<-ERCA,HR -3136.739214776165\n",
       "HRSA<-ERCA,HR -3236.139444268571\n",
       "HYP<-APL -10075.035674773593\n",
       "INT<- -6658.474123732305\n",
       "KINK<-VLNG,VTUB -2857.4047588872545\n",
       "VTUB<-DISC,VMCH -3418.0624644648706\n",
       "MINV<-INT,VLNG -3909.141376477159\n",
       "MVS<- -8115.930179708217\n",
       "PAP<- -8357.309503322173\n",
       "PCWP<-LVV -4417.152413059788\n",
       "PMB<-PAP -782.1722102996109\n",
       "PRSS<-INT,KINK,VTUB -17104.845843232193\n",
       "PVS<-FIO2,VALV -1881.4613201566744\n",
       "SHNT<-INT,PMB -3967.786877754141\n",
       "VMCH<-MVS -6171.55480017248\n",
       "**********\n",
       "bnlearn modelstring = \n",
       "[ACO2|VALV:CCHL][VALV|INT:VLNG][CCHL|TPR:SAO2][ANES][APL][BP|TPR:CO][TPR|APL][CO|HR:STKV][SAO2|SHNT:PVS][HR|CCHL][STKV|LVF:HYP][CVP|LVV][LVV|LVF:HYP][DISC][ECO2|VLNG:ACO2][VLNG|INT:VTUB][ERCA][ERLO][FIO2][HIST|LVF][LVF][HRBP|ERLO:HR][HREK|ERCA:HR][HRSA|ERCA:HR][HYP|APL][INT][KINK|VLNG:VTUB][VTUB|DISC:VMCH][MINV|INT:VLNG][MVS][PAP][PCWP|LVV][PMB|PAP][PRSS|INT:VTUB:KINK][PVS|VALV:FIO2][SHNT|INT:PMB][VMCH|MVS]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "library(reticulate)\n",
    "m <- import(\"pygobnilp.gobnilp\")$Gobnilp()\n",
    "m$learn(alarm,plot=FALSE,score='DiscreteBIC')\n",
    "m$learned_bn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gobnilp is a lot slower, but on the other hand, the learned BN has a much better BIC score of -218632.30999687716. Note that Gobnilp is being run here with its default limit of at most 3 parents for any node. It could be that an even higher-scoring network exists where some nodes have more than 3 parents.\n",
    "\n",
    "We can get bnlearn to score the network learned by Gobnilp to check that we get the same BIC score."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "-218632.315679648"
      ],
      "text/latex": [
       "-218632.315679648"
      ],
      "text/markdown": [
       "-218632.315679648"
      ],
      "text/plain": [
       "[1] -218632.3"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "score(model2network(m$learned_bn$bnlearn_modelstring()),alarm,type=\"bic\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sure enough we get the same score (modulo some numerical imprecision). Now let's run the same experiment with bnlearn's built-in hailfinder dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "  Bayesian network learned via Score-based methods\n",
       "\n",
       "  model:\n",
       "   [N07muVerMo][SubjVertMo][QGVertMotion][SatContMoist][RaoContMoist]\n",
       "   [VISCloudCov][IRCloudCover][AMInstabMt][WndHodograph][MorningBound]\n",
       "   [LoLevMoistAd][Date][MorningCIN][LIfr12ZDENSd][AMDewptCalPl][LatestCIN][LLIW]\n",
       "   [CombVerMo|N07muVerMo:SubjVertMo:QGVertMotion]\n",
       "   [CombMoisture|SatContMoist:RaoContMoist][CombClouds|VISCloudCov:IRCloudCover]\n",
       "   [Scenario|Date][CurPropConv|LatestCIN:LLIW][AreaMesoALS|CombVerMo]\n",
       "   [AreaMoDryAir|CombVerMo:CombMoisture][ScenRelAMCIN|Scenario]\n",
       "   [ScenRelAMIns|Scenario][ScenRel34|Scenario][ScnRelPlFcst|Scenario]\n",
       "   [Dewpoints|Scenario][LowLLapse|Scenario][MeanRH|Scenario][MidLLapse|Scenario]\n",
       "   [MvmtFeatures|Scenario][RHRatio|Scenario][SfcWndShfDis|Scenario]\n",
       "   [SynForcng|Scenario][TempDis|Scenario][WindAloft|Scenario]\n",
       "   [WindFieldMt|Scenario][WindFieldPln|Scenario]\n",
       "   [CldShadeOth|CombVerMo:AreaMoDryAir:CombClouds]\n",
       "   [AMCINInScen|ScenRelAMCIN:MorningCIN]\n",
       "   [AMInsWliScen|ScenRelAMIns:LIfr12ZDENSd:AMDewptCalPl]\n",
       "   [InsInMt|CldShadeOth:AMInstabMt][OutflowFrMt|InsInMt:WndHodograph]\n",
       "   [CldShadeConv|InsInMt:WndHodograph][MountainFcst|InsInMt]\n",
       "   [Boundaries|WndHodograph:OutflowFrMt:MorningBound]\n",
       "   [CompPlFcst|CombVerMo:CldShadeOth:CldShadeConv][CapChange|CompPlFcst]\n",
       "   [InsChange|CompPlFcst:LoLevMoistAd][CapInScen|CompPlFcst:AMCINInScen]\n",
       "   [InsSclInScen|InsChange:AMInsWliScen]\n",
       "   [PlainsFcst|ScenRelAMCIN:InsSclInScen:CurPropConv]\n",
       "   [N34StarFcst|ScenRel34:PlainsFcst][R5Fcst|MountainFcst:N34StarFcst]\n",
       "  nodes:                                 56 \n",
       "  arcs:                                  64 \n",
       "    undirected arcs:                     0 \n",
       "    directed arcs:                       64 \n",
       "  average markov blanket size:           3.25 \n",
       "  average neighbourhood size:            2.29 \n",
       "  average branching factor:              1.14 \n",
       "\n",
       "  learning algorithm:                    Hill-Climbing \n",
       "  score:                                 BIC (disc.) \n",
       "  penalization coefficient:              4.951744 \n",
       "  tests used in the learning procedure:  5060 \n",
       "  optimized:                             TRUE \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "-990474.753096428"
      ],
      "text/latex": [
       "-990474.753096428"
      ],
      "text/markdown": [
       "-990474.753096428"
      ],
      "text/plain": [
       "[1] -990474.8"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "  Bayesian network learned via Score-based methods\n",
       "\n",
       "  model:\n",
       "   [N07muVerMo][SubjVertMo][QGVertMotion][SatContMoist][RaoContMoist]\n",
       "   [VISCloudCov][IRCloudCover][AMInstabMt][WndHodograph][MorningBound]\n",
       "   [LoLevMoistAd][Date][MorningCIN][LIfr12ZDENSd][AMDewptCalPl][LatestCIN][LLIW]\n",
       "   [CombVerMo|N07muVerMo:SubjVertMo:QGVertMotion]\n",
       "   [CombMoisture|SatContMoist:RaoContMoist][CombClouds|VISCloudCov:IRCloudCover]\n",
       "   [Scenario|Date][CurPropConv|LatestCIN:LLIW][AreaMesoALS|CombVerMo]\n",
       "   [AreaMoDryAir|CombVerMo:CombMoisture][ScenRelAMCIN|Scenario]\n",
       "   [ScenRelAMIns|Scenario][ScenRel34|Scenario][ScnRelPlFcst|Scenario]\n",
       "   [Dewpoints|Scenario][LowLLapse|Scenario][MeanRH|Scenario][MidLLapse|Scenario]\n",
       "   [MvmtFeatures|Scenario][RHRatio|Scenario][SfcWndShfDis|Scenario]\n",
       "   [SynForcng|Scenario][TempDis|Scenario][WindAloft|Scenario]\n",
       "   [WindFieldMt|Scenario][WindFieldPln|Scenario]\n",
       "   [CldShadeOth|CombVerMo:AreaMoDryAir:CombClouds]\n",
       "   [AMCINInScen|ScenRelAMCIN:MorningCIN]\n",
       "   [AMInsWliScen|ScenRelAMIns:LIfr12ZDENSd:AMDewptCalPl]\n",
       "   [InsInMt|CldShadeOth:AMInstabMt][OutflowFrMt|InsInMt:WndHodograph]\n",
       "   [CldShadeConv|InsInMt:WndHodograph][MountainFcst|InsInMt]\n",
       "   [Boundaries|WndHodograph:OutflowFrMt:MorningBound]\n",
       "   [CompPlFcst|CombVerMo:CldShadeOth:CldShadeConv][CapChange|CompPlFcst]\n",
       "   [InsChange|CompPlFcst:LoLevMoistAd][CapInScen|CompPlFcst:AMCINInScen]\n",
       "   [InsSclInScen|InsChange:AMInsWliScen]\n",
       "   [PlainsFcst|ScenRelAMCIN:InsSclInScen:CurPropConv]\n",
       "   [N34StarFcst|ScenRel34:PlainsFcst][R5Fcst|MountainFcst:N34StarFcst]\n",
       "  nodes:                                 56 \n",
       "  arcs:                                  64 \n",
       "    undirected arcs:                     0 \n",
       "    directed arcs:                       64 \n",
       "  average markov blanket size:           3.25 \n",
       "  average neighbourhood size:            2.29 \n",
       "  average branching factor:              1.14 \n",
       "\n",
       "  learning algorithm:                    Tabu Search \n",
       "  score:                                 BIC (disc.) \n",
       "  penalization coefficient:              4.951744 \n",
       "  tests used in the learning procedure:  6014 \n",
       "  optimized:                             TRUE \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "-990474.753096428"
      ],
      "text/latex": [
       "-990474.753096428"
      ],
      "text/markdown": [
       "-990474.753096428"
      ],
      "text/plain": [
       "[1] -990474.8"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "**********\n",
       "BN has score -990162.8050400699\n",
       "**********\n",
       "AMCINInScen<-MorningCIN,ScenRelAMCIN -14827.804016810527\n",
       "ScenRelAMCIN<-Scenario -54.46918153893982\n",
       "MorningCIN<- -22451.66219109244\n",
       "AMDewptCalPl<- -21401.356146585407\n",
       "AMInsWliScen<-AMDewptCalPl,ScenRelAMIns -19662.46130603437\n",
       "ScenRelAMIns<-ScnRelPlFcst -272.34590769473465\n",
       "AMInstabMt<- -21979.829675807607\n",
       "AreaMesoALS<-CombVerMo -59.42092531521676\n",
       "CombVerMo<-N07muVerMo,QGVertMotion,SubjVertMo -16140.95428595368\n",
       "AreaMoDryAir<-CombMoisture,CombVerMo -15243.049284360612\n",
       "CombMoisture<-RaoContMoist,SatContMoist -22201.123014157503\n",
       "Boundaries<-MorningBound,OutflowFrMt,WndHodograph -15464.935101862531\n",
       "OutflowFrMt<-InsInMt,WndHodograph -12684.806287718633\n",
       "WndHodograph<- -27534.501610060404\n",
       "MorningBound<- -19880.990986381013\n",
       "CapChange<-CompPlFcst -29.71046265760838\n",
       "CompPlFcst<-CldShadeConv,CldShadeOth,CombVerMo -20195.300816671996\n",
       "CapInScen<-AMCINInScen,CapChange -8769.900372311302\n",
       "CldShadeConv<-InsInMt,WndHodograph -11930.685819990145\n",
       "InsInMt<-AMInstabMt,CldShadeOth -9186.457382838727\n",
       "CldShadeOth<-AreaMoDryAir,CombClouds,CombVerMo -10808.429735453852\n",
       "CombClouds<-IRCloudCover,VISCloudCov -12416.50118390726\n",
       "VISCloudCov<- -18807.3986119801\n",
       "IRCloudCover<- -20081.087288850224\n",
       "RaoContMoist<- -26451.36187705046\n",
       "SatContMoist<- -26381.105022929798\n",
       "QGVertMotion<- -24641.079933873705\n",
       "SubjVertMo<- -24741.633744298222\n",
       "N07muVerMo<- -27740.002999058954\n",
       "CurPropConv<-PlainsFcst -26471.80989062072\n",
       "PlainsFcst<-InsSclInScen,Scenario -17466.172736436805\n",
       "Date<-Scenario -33427.328674309276\n",
       "Scenario<-ScnRelPlFcst -544.6918153894782\n",
       "Dewpoints<-Scenario -30495.957816711856\n",
       "InsChange<-CompPlFcst,LoLevMoistAd -15438.63271567674\n",
       "LoLevMoistAd<- -26697.017932815765\n",
       "InsSclInScen<-AMInsWliScen,InsChange -11294.899810438172\n",
       "LIfr12ZDENSd<-AMDewptCalPl,AMInsWliScen -21266.227971706914\n",
       "LLIW<-CurPropConv -21747.034006031623\n",
       "LatestCIN<-CurPropConv,LLIW -19532.950433050428\n",
       "LowLLapse<-ScnRelPlFcst -22799.692963671845\n",
       "ScnRelPlFcst<-WindFieldMt -43976.15779083476\n",
       "MeanRH<-Scenario -16858.5681959753\n",
       "MidLLapse<-Scenario -18464.904348591688\n",
       "MountainFcst<-InsInMt -15270.69659773143\n",
       "MvmtFeatures<-Scenario -19160.848763138114\n",
       "N34StarFcst<-PlainsFcst,ScenRel34 -5277.188589205412\n",
       "ScenRel34<-ScnRelPlFcst -217.8767261557948\n",
       "R5Fcst<-MountainFcst,N34StarFcst -89.1313879728207\n",
       "RHRatio<-ScnRelPlFcst -16934.816498959048\n",
       "WindFieldMt<- -13845.133818871916\n",
       "SfcWndShfDis<-ScnRelPlFcst -28237.61080138319\n",
       "SynForcng<-ScnRelPlFcst -27008.48382226494\n",
       "TempDis<-Scenario -20415.645867981773\n",
       "WindAloft<-ScnRelPlFcst -19892.419524674344\n",
       "WindFieldPln<-Scenario -25290.54036622376\n",
       "**********\n",
       "bnlearn modelstring = \n",
       "[AMCINInScen|ScenRelAMCIN:MorningCIN][ScenRelAMCIN|Scenario][MorningCIN][AMDewptCalPl][AMInsWliScen|ScenRelAMIns:AMDewptCalPl][ScenRelAMIns|ScnRelPlFcst][AMInstabMt][AreaMesoALS|CombVerMo][CombVerMo|QGVertMotion:SubjVertMo:N07muVerMo][AreaMoDryAir|CombVerMo:CombMoisture][CombMoisture|RaoContMoist:SatContMoist][Boundaries|OutflowFrMt:WndHodograph:MorningBound][OutflowFrMt|WndHodograph:InsInMt][WndHodograph][MorningBound][CapChange|CompPlFcst][CompPlFcst|CldShadeConv:CombVerMo:CldShadeOth][CapInScen|CapChange:AMCINInScen][CldShadeConv|WndHodograph:InsInMt][InsInMt|AMInstabMt:CldShadeOth][CldShadeOth|AreaMoDryAir:CombClouds:CombVerMo][CombClouds|VISCloudCov:IRCloudCover][VISCloudCov][IRCloudCover][RaoContMoist][SatContMoist][QGVertMotion][SubjVertMo][N07muVerMo][CurPropConv|PlainsFcst][PlainsFcst|Scenario:InsSclInScen][Date|Scenario][Scenario|ScnRelPlFcst][Dewpoints|Scenario][InsChange|LoLevMoistAd:CompPlFcst][LoLevMoistAd][InsSclInScen|InsChange:AMInsWliScen][LIfr12ZDENSd|AMDewptCalPl:AMInsWliScen][LLIW|CurPropConv][LatestCIN|CurPropConv:LLIW][LowLLapse|ScnRelPlFcst][ScnRelPlFcst|WindFieldMt][MeanRH|Scenario][MidLLapse|Scenario][MountainFcst|InsInMt][MvmtFeatures|Scenario][N34StarFcst|ScenRel34:PlainsFcst][ScenRel34|ScnRelPlFcst][R5Fcst|N34StarFcst:MountainFcst][RHRatio|ScnRelPlFcst][WindFieldMt][SfcWndShfDis|ScnRelPlFcst][SynForcng|ScnRelPlFcst][TempDis|Scenario][WindAloft|ScnRelPlFcst][WindFieldPln|Scenario]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "-990162.80504007"
      ],
      "text/latex": [
       "-990162.80504007"
      ],
      "text/markdown": [
       "-990162.80504007"
      ],
      "text/plain": [
       "[1] -990162.8"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data(hailfinder)\n",
    "hailfinderhc.bn <- hc(hailfinder)\n",
    "hailfinderhc.bn\n",
    "score(hailfinderhc.bn,hailfinder,type=\"bic\")\n",
    "hailfindertabu.bn <- tabu(hailfinder)\n",
    "hailfindertabu.bn\n",
    "score(hailfindertabu.bn,hailfinder,type=\"bic\")\n",
    "m$learn(hailfinder,plot=FALSE,score='DiscreteBIC')\n",
    "m$learned_bn\n",
    "score(model2network(m$learned_bn$bnlearn_modelstring()),hailfinder,type=\"bic\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gobnilp takes a very long time to find a BN mainly due to this Python version of Gobnilp computing the necessary local scores slowly. The found BN has a BIC score of -990162.80504007 in contrast to -990474.753096428 BIC score achieved by hill-climbing and tabu.\n",
    "\n",
    "Finally, let's have a look at the Insurance dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "  Bayesian network learned via Score-based methods\n",
       "\n",
       "  model:\n",
       "   [RuggedAuto][MakeModel|RuggedAuto][CarValue|RuggedAuto:MakeModel]\n",
       "   [Mileage|CarValue][VehicleYear|MakeModel:Mileage:CarValue]\n",
       "   [SocioEcon|VehicleYear:MakeModel][Antilock|VehicleYear:MakeModel]\n",
       "   [Airbag|VehicleYear:MakeModel][ThisCarDam|Mileage:Antilock]\n",
       "   [OtherCar|SocioEcon][Cushioning|RuggedAuto:Airbag]\n",
       "   [Accident|ThisCarDam:RuggedAuto][ThisCarCost|ThisCarDam:CarValue]\n",
       "   [MedCost|ThisCarDam:Cushioning][DrivQuality|Accident:Mileage]\n",
       "   [Theft|ThisCarDam:ThisCarCost][OtherCarCost|RuggedAuto:Accident]\n",
       "   [ILiCost|Accident][Age|SocioEcon:DrivQuality]\n",
       "   [PropCost|ThisCarCost:OtherCarCost][GoodStudent|Age:SocioEcon]\n",
       "   [SeniorTrain|Age:DrivQuality][RiskAversion|Age:DrivQuality:SeniorTrain]\n",
       "   [DrivingSkill|RiskAversion:DrivQuality][HomeBase|SocioEcon:RiskAversion]\n",
       "   [AntiTheft|SocioEcon:RiskAversion][DrivHist|RiskAversion:DrivingSkill]\n",
       "  nodes:                                 27 \n",
       "  arcs:                                  50 \n",
       "    undirected arcs:                     0 \n",
       "    directed arcs:                       50 \n",
       "  average markov blanket size:           4.44 \n",
       "  average neighbourhood size:            3.70 \n",
       "  average branching factor:              1.85 \n",
       "\n",
       "  learning algorithm:                    Hill-Climbing \n",
       "  score:                                 BIC (disc.) \n",
       "  penalization coefficient:              4.951744 \n",
       "  tests used in the learning procedure:  1755 \n",
       "  optimized:                             TRUE \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "-266113.028472967"
      ],
      "text/latex": [
       "-266113.028472967"
      ],
      "text/markdown": [
       "-266113.028472967"
      ],
      "text/plain": [
       "[1] -266113"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "\n",
       "  Bayesian network learned via Score-based methods\n",
       "\n",
       "  model:\n",
       "   [MakeModel][CarValue|MakeModel][Mileage|CarValue]\n",
       "   [VehicleYear|MakeModel:Mileage:CarValue][SocioEcon|VehicleYear:MakeModel]\n",
       "   [RuggedAuto|VehicleYear:MakeModel][Antilock|VehicleYear:MakeModel]\n",
       "   [Airbag|VehicleYear:MakeModel][ThisCarDam|Mileage:Antilock]\n",
       "   [OtherCar|SocioEcon][Cushioning|RuggedAuto:Airbag]\n",
       "   [Accident|ThisCarDam:RuggedAuto][ThisCarCost|ThisCarDam:CarValue]\n",
       "   [MedCost|ThisCarDam:Cushioning][DrivQuality|Accident:Mileage]\n",
       "   [Theft|ThisCarDam:ThisCarCost][OtherCarCost|RuggedAuto:Accident]\n",
       "   [ILiCost|Accident][Age|SocioEcon:DrivQuality]\n",
       "   [PropCost|ThisCarCost:OtherCarCost][GoodStudent|Age:SocioEcon]\n",
       "   [SeniorTrain|Age:DrivQuality][RiskAversion|Age:DrivQuality:SeniorTrain]\n",
       "   [DrivingSkill|RiskAversion:DrivQuality][HomeBase|SocioEcon:RiskAversion]\n",
       "   [AntiTheft|SocioEcon:RiskAversion][DrivHist|RiskAversion:DrivingSkill]\n",
       "  nodes:                                 27 \n",
       "  arcs:                                  50 \n",
       "    undirected arcs:                     0 \n",
       "    directed arcs:                       50 \n",
       "  average markov blanket size:           4.44 \n",
       "  average neighbourhood size:            3.70 \n",
       "  average branching factor:              1.85 \n",
       "\n",
       "  learning algorithm:                    Tabu Search \n",
       "  score:                                 BIC (disc.) \n",
       "  penalization coefficient:              4.951744 \n",
       "  tests used in the learning procedure:  2405 \n",
       "  optimized:                             TRUE \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "-265500.923940094"
      ],
      "text/latex": [
       "-265500.923940094"
      ],
      "text/markdown": [
       "-265500.923940094"
      ],
      "text/plain": [
       "[1] -265500.9"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "**********\n",
       "BN has score -264708.85975077713\n",
       "**********\n",
       "Accident<-Antilock,DrivQuality -12310.195471829667\n",
       "DrivQuality<-DrivingSkill,RiskAversion -4890.4371039384205\n",
       "Antilock<-MakeModel,VehicleYear -3057.5537899833384\n",
       "Age<-RiskAversion,SocioEcon -17309.67690764885\n",
       "RiskAversion<- -22128.556172909437\n",
       "SocioEcon<-RiskAversion -21914.967878236253\n",
       "Airbag<-MakeModel,VehicleYear -4426.614280155772\n",
       "VehicleYear<-SocioEcon -10055.937028683162\n",
       "MakeModel<-SocioEcon -18409.07715172664\n",
       "AntiTheft<-RiskAversion,SocioEcon -5427.903744351533\n",
       "CarValue<-MakeModel,Mileage,VehicleYear -11811.232170180529\n",
       "Mileage<-Accident,DrivQuality -23647.681978155444\n",
       "Cushioning<-Airbag,RuggedAuto -14785.848505297057\n",
       "RuggedAuto<-MakeModel,VehicleYear -12165.511166709075\n",
       "DrivHist<-DrivingSkill,RiskAversion -9560.920896777798\n",
       "DrivingSkill<-Age,SeniorTrain -17818.740523142405\n",
       "SeniorTrain<-Age,RiskAversion -1712.9733331992413\n",
       "GoodStudent<-Age,SocioEcon -1937.6149662237774\n",
       "HomeBase<-RiskAversion,SocioEcon -14317.997554087879\n",
       "ILiCost<-Accident -2459.8973674972813\n",
       "MedCost<-Accident,Cushioning -3741.377927855121\n",
       "OtherCar<-SocioEcon -10846.63435568069\n",
       "OtherCarCost<-Accident,RuggedAuto -4190.349482680152\n",
       "PropCost<-OtherCarCost,ThisCarCost -11581.005875670833\n",
       "ThisCarCost<-CarValue,ThisCarDam -1499.5192944884964\n",
       "Theft<-ThisCarCost,ThisCarDam -120.29650307184178\n",
       "ThisCarDam<-Accident,RuggedAuto -2580.3383205964406\n",
       "**********\n",
       "bnlearn modelstring = \n",
       "[Accident|DrivQuality:Antilock][DrivQuality|RiskAversion:DrivingSkill][Antilock|VehicleYear:MakeModel][Age|RiskAversion:SocioEcon][RiskAversion][SocioEcon|RiskAversion][Airbag|VehicleYear:MakeModel][VehicleYear|SocioEcon][MakeModel|SocioEcon][AntiTheft|RiskAversion:SocioEcon][CarValue|VehicleYear:MakeModel:Mileage][Mileage|DrivQuality:Accident][Cushioning|RuggedAuto:Airbag][RuggedAuto|VehicleYear:MakeModel][DrivHist|RiskAversion:DrivingSkill][DrivingSkill|SeniorTrain:Age][SeniorTrain|RiskAversion:Age][GoodStudent|Age:SocioEcon][HomeBase|RiskAversion:SocioEcon][ILiCost|Accident][MedCost|Cushioning:Accident][OtherCar|SocioEcon][OtherCarCost|Accident:RuggedAuto][PropCost|OtherCarCost:ThisCarCost][ThisCarCost|ThisCarDam:CarValue][Theft|ThisCarCost:ThisCarDam][ThisCarDam|Accident:RuggedAuto]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "-264708.859750777"
      ],
      "text/latex": [
       "-264708.859750777"
      ],
      "text/markdown": [
       "-264708.859750777"
      ],
      "text/plain": [
       "[1] -264708.9"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "data(insurance)\n",
    "insurancehc.bn <- hc(insurance)\n",
    "insurancehc.bn\n",
    "score(insurancehc.bn,insurance,type=\"bic\")\n",
    "insurancetabu.bn <- tabu(insurance)\n",
    "insurancetabu.bn\n",
    "score(insurancetabu.bn,insurance,type=\"bic\")\n",
    "m$learn(insurance,plot=FALSE,score='DiscreteBIC')\n",
    "m$learned_bn\n",
    "score(model2network(m$learned_bn$bnlearn_modelstring()),insurance,type=\"bic\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So hill-climbing finds a BN with BIC score -266113.028472967, tabu search manages -265500.923940094 and Gobnilp achieves -264708.85975077713. In this case Gobnilp was not too slow either (since Insurance only has 27 nodes and the default limit of at most 3 parents per node is in use.)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.4.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
